////////////////////////////////////////////////////////////////////// 
// libsrc/PedigreeGlobals.cpp 
// (c) 2000-2008 Goncalo Abecasis
// 
// This file is distributed as part of the MaCH source code package   
// and may not be redistributed in any form, without prior written    
// permission from the author. Permission is granted for you to       
// modify this file for your own personal use, but modified versions  
// must retain this copyright notice and must not be distributed.     
// 
// Permission is granted for you to use this file to compile MaCH.    
// 
// All computer programs have bugs. Use this file at your own risk.   
// 
// Saturday April 12, 2008
// 
 
#include "PedigreeGlobals.h"
#include "Sort.h"
#include "Error.h"

#include <math.h>
#include <string.h>
#include <ctype.h>

int PedigreeGlobals::traitCount = 0;
int PedigreeGlobals::affectionCount = 0;
int PedigreeGlobals::covariateCount = 0;
int PedigreeGlobals::markerCount = 0;

// If this value isn't set, all X chromosome data will be rejected
bool PedigreeGlobals::chromosomeX = false;
bool PedigreeGlobals::sexSpecificMap = false;

StringArray   PedigreeGlobals::traitNames;
StringArray   PedigreeGlobals::markerNames;
StringArray   PedigreeGlobals::covariateNames;
StringArray   PedigreeGlobals::affectionNames;
StringIntHash PedigreeGlobals::markerLookup;
StringIntHash PedigreeGlobals::traitLookup;
StringIntHash PedigreeGlobals::affectionLookup;
StringIntHash PedigreeGlobals::covariateLookup;

int PedigreeGlobals::markerInfoCount = 0;
int PedigreeGlobals::markerInfoSize = 0;

MarkerInfo ** PedigreeGlobals::markerInfo = NULL;
MarkerInfo ** PedigreeGlobals::markerInfoByInteger = NULL;
StringHash    PedigreeGlobals::markerInfoByName;

int MarkerInfo::count = 0;

int MarkerInfo::ComparePosition(MarkerInfo ** left, MarkerInfo ** right)
   {
   if ((*left)->chromosome != (*right)->chromosome)
      return (*left)->chromosome - (*right)->chromosome;

   double difference = (*left)->position - (*right)->position;

   if (difference >  0.0)
      return 1;
   else if (difference == 0.0)
      return (*left)->serial - (*right)->serial;
   else
      return -1;
   }

String MarkerInfo::GetAlleleLabel(int allele) 
   {
   if (alleleLabels.Length() > allele && alleleLabels[allele].Length())
      return alleleLabels[allele];
   else if (alleleLabels.Length() <= allele)
      alleleLabels.Dimension(allele + 1);
   return alleleLabels[allele] = allele;
   }

bool MarkerInfo::AdjustFrequencies()
   {
   if (freq.Length() <= 1)
      {
      freq.Clear();
      return false;
      }

   if (freq.Min() < 0.0)
      error("Locus %s has negative allele frequencies\n", (const char *) name);

   double sum = freq.Sum();

   if (sum <= 0.0)
      error("Locus %s frequencies sum to %f, which doesn't make sense\n",
         (const char *) name, sum);

   if ( sum != 1.0)
      freq *= 1.0 / sum;

   if ( fabs(sum - 1.0) > 1.2e-5 )
      {
      printf("Locus %s frequencies sum to %f, adjusted to 1.0\n",
          (const char *) name, sum);
      return true;
      }

   return false;
   }

void MarkerInfo::IndexAlleles()
   {
   if (alleleLabels.Length() >= 255)
      error("Marker %s has more than 254 distinct alleles\n",
            (const char *) name);

   alleleNumbers.Clear();

   for (int i = 1; i < alleleLabels.Length(); i++)
      alleleNumbers.SetInteger(alleleLabels[i], i);
   }

int MarkerInfo::NewAllele(const String & label)
   {
   if (alleleLabels.Length() == 0)
      alleleLabels.Push("");

   if (alleleLabels.Length() >= 255)
      error("Marker %s has more than 254 distinct alleles\n",
            (const char *) name);

   alleleNumbers.SetInteger(label, alleleLabels.Length());
   alleleLabels.Push(label);

   return alleleLabels.Length() - 1;
   }

int PedigreeGlobals::GetTraitID(const char * name)
   {
   int idx = traitLookup.Integer(name);

   if (idx != -1) return idx;

   traitNames.Add(name);
   traitLookup.SetInteger(name, traitCount);
   return traitCount++;
   }

int PedigreeGlobals::GetAffectionID(const char * name)
   {
   int idx = affectionLookup.Integer(name);

   if (idx != -1) return idx;

   affectionNames.Add(name);
   affectionLookup.SetInteger(name, affectionCount);
   return affectionCount++;
   }

int PedigreeGlobals::GetCovariateID(const char * name)
   {
   int idx = covariateLookup.Integer(name);

   if (idx != -1) return idx;

   covariateNames.Add(name);
   covariateLookup.SetInteger(name, covariateCount);
   return covariateCount++;
   }

int PedigreeGlobals::GetMarkerID(const char * name)
   {
   int idx = markerLookup.Integer(name);

   if (idx != -1) return idx;

   markerNames.Add(name);
   markerLookup.SetInteger(name, markerCount);

   // Grow the marker info key ...
   if (markerCount == 0)
      {
      markerInfoByInteger = new MarkerInfo * [16];

      for (int i = 0; i < 16; i++)
         markerInfoByInteger[i] = NULL;
      }
   else if ((markerCount & (markerCount - 1)) == 0 && markerCount > 15)
      {
      MarkerInfo ** newKey = new MarkerInfo * [markerCount * 2];

      for (int i = 0; i < markerCount; i++)
         newKey[i] = markerInfoByInteger[i];

      for (int i = markerCount; i < markerCount * 2; i++)
         newKey[i] = NULL;

      delete [] markerInfoByInteger;

      markerInfoByInteger = newKey;
      }

   return markerCount++;
   }

MarkerInfo * PedigreeGlobals::GetMarkerInfo(String & name)
   {
   MarkerInfo * info = (MarkerInfo *) markerInfoByName.Object(name);

   if (info != NULL) return info;

   info = new MarkerInfo(name);
   markerInfoByName.Add(name, info);

   if (markerInfoCount >= markerInfoSize)
      GrowMarkerInfo();

   markerInfo[markerInfoCount++] = info;

   int markerId = LookupMarker(name);
   if (markerId >= 0) markerInfoByInteger[markerId] = info;

   return info;
   }

MarkerInfo * PedigreeGlobals::GetMarkerInfo(int markerId)
   {
   if (markerId >= markerCount)
      error("Attempted to retrieve MarkerInfo using out-of-bounds index\n");

   if (markerInfoByInteger[markerId] != NULL)
      return markerInfoByInteger[markerId];
   else
      return GetMarkerInfo(markerNames[markerId]);
   }

void PedigreeGlobals::GrowMarkerInfo()
   {
   int newSize = markerInfoSize ? 2 * markerInfoSize : 32;

   MarkerInfo ** newArray = new MarkerInfo * [newSize];

   if (markerInfoSize)
      {
      memcpy(newArray, markerInfo, sizeof(MarkerInfo *) * markerInfoSize);
      delete [] markerInfo;
      }

   markerInfo = newArray;
   markerInfoSize = newSize;
   }

void PedigreeGlobals::FlagMissingMarkers(IntArray & missingMarkers)
   {
   int skipped_markers = 0;

   if (missingMarkers.Length())
      {
      StringArray names;

      printf("These markers couldn't be placed and won't be analysed:");

      for (int i = 0; i < missingMarkers.Length(); i++)
         names.Push(GetMarkerInfo(missingMarkers[i])->name);
      names.Sort();

      for (int i = 0, line = 80, lines = 0; i < missingMarkers.Length(); i++)
         {
         if (line + names[i].Length() + 1 > 79)
            printf("\n   "), line = 3, lines++;

         if (lines < 5)
            {
            printf("%s ", (const char *) names[i]);
            line += names[i].Length() + 1;
            }
         else
            skipped_markers++;
         }

      if (skipped_markers)
         printf("as well as %d other unlisted markers...", skipped_markers);

      printf("\n\n");
      }
   }

void PedigreeGlobals::GetOrderedMarkers(IntArray & markers)
   {
   if (markers.Length() == 0)
      {
      markers.Dimension(markerCount);
      markers.SetSequence(0, 1);
      }

   MarkerInfo ** subset = new MarkerInfo * [markers.Length()];

   int count = 0;
   IntArray missingMarkers;

   for (int i = 0; i < markers.Length(); i++)
      {
      MarkerInfo * info = GetMarkerInfo(markers[i]);

      if (info->chromosome != -1)
         subset[count++] = info;
      else
         missingMarkers.Push(i);
      }

   FlagMissingMarkers(missingMarkers);

   QuickSort(subset, count, sizeof(MarkerInfo *),
             COMPAREFUNC MarkerInfo::ComparePosition);

   markers.Clear();
   for (int i = 0; i < count; i++)
      markers.Push(GetMarkerID(subset[i]->name));
   }

int PedigreeGlobals::SortMarkersInMapOrder(IntArray & markers, int chromosome)
   {
   if (markers.Length() == 0)
      {
      markers.Dimension(markerCount);
      markers.SetSequence(0, 1);
      }

   MarkerInfo ** subset = new MarkerInfo * [markers.Length()];

   int count = 0;
   IntArray missingMarkers;

   for (int i = 0; i < markers.Length(); i++)
      {
      MarkerInfo * info = GetMarkerInfo(markers[i]);

      if (info->chromosome != -1)
         subset[count++] = info;
      else if (chromosome == -1)
         missingMarkers.Push(i);
      }

   if (chromosome == -1)
      FlagMissingMarkers(missingMarkers);

   QuickSort(subset, count, sizeof(MarkerInfo *),
             COMPAREFUNC MarkerInfo::ComparePosition);

   markers.Clear();

   int  current_chromosome = -1, next_chromosome = 0;

   for (int i = 0; i < count; i++)
      if (subset[i]->chromosome < chromosome)
         continue;
      else if (current_chromosome == -1 ||
               subset[i]->chromosome == current_chromosome)
         {
         markers.Push(GetMarkerID(subset[i]->name));
         current_chromosome = subset[i]->chromosome;
         }
      else if (!next_chromosome)
         {
         next_chromosome = subset[i]->chromosome;
         break;
         }

   delete [] subset;

   return next_chromosome;
   }

void PedigreeGlobals::VerifySexSpecificOrder()
   {
   if (markerCount <= 1)
      return;

   MarkerInfo ** sortedMarkers = new MarkerInfo * [markerCount];

   for (int i = 0; i < markerCount; i++)
      sortedMarkers[i] = GetMarkerInfo(i);

   QuickSort(sortedMarkers, markerCount, sizeof(MarkerInfo *),
             COMPAREFUNC MarkerInfo::ComparePosition);

   double prev_female = sortedMarkers[0]->positionFemale;
   double prev_male = sortedMarkers[0]->positionMale;
   double curr_female, curr_male;

   int    prev_chromosome = sortedMarkers[0]->chromosome;
   int    curr_chromosome;

   for (int i = 1; i < markerCount; i++)
      {
      curr_chromosome = sortedMarkers[i]->chromosome;
      curr_female = sortedMarkers[i]->positionFemale;
      curr_male = sortedMarkers[i]->positionMale;

      if (curr_chromosome == prev_chromosome &&
         (curr_female < prev_female || curr_male < prev_male))
         error("Sex-specific and sex-averaged maps are inconsistent.\n\n"
               "In the sex-averaged map, marker %s (%.2f cM) follows marker %s (%.2f cM).\n"
               "In the %smale map, marker %s (%.2f cM) PRECEDES marker %s (%.2f cM).\n",
               (const char *) sortedMarkers[i]->name,
               sortedMarkers[i]->position * 100,
               (const char *) sortedMarkers[i-1]->name,
               sortedMarkers[i-1]->position * 100,
               curr_female < prev_female ? "fe" : "",
               (const char *) sortedMarkers[i]->name,
               (curr_female < prev_female ? curr_female : curr_male) * 100,
               (const char *) sortedMarkers[i-1]->name,
               (curr_female < prev_female ? prev_female : prev_male) * 100);

      prev_chromosome = curr_chromosome;
      prev_female = curr_female;
      prev_male = curr_male;
      }

   delete [] sortedMarkers;
   }

void PedigreeGlobals::LoadAlleleFrequencies(const char * filename, bool required)
   {
   // This function is often called with an empty string, and not
   // all implementations of the C library like that ...
   if (filename[0] == 0)
      if (required)
         error("No name provided for required allele freuquency file\n");
      else
         return;

   // If we get here, the filename is not empty and things should
   // work as planned
   IFILE f = ifopen(filename, "rb");

   if (f == NULL)
      if (required)
         error("Failed to open required alallele frequency file '%s'",
               (const char *) filename);
      else
         return;

   LoadAlleleFrequencies(f);
   ifclose(f);
   }

void PedigreeGlobals::LoadAlleleFrequencies(IFILE & input)
   {
   int         done = 0;
   String      buffer;
   StringArray tokens;
   MarkerInfo *info = NULL;

   bool need_blank_line = false;
   int  allele_size, old_max, next_allele = 0; // Initialization avoids compiler warning

   while (!ifeof(input) && !done)
      {
      int   i, j;

      buffer.ReadLine(input);

      tokens.Clear();
      tokens.AddTokens(buffer, WHITESPACE);

      if (tokens.Length() < 1) continue;

      switch (toupper(tokens[0][0]))
         {
         case 'M' :
            if (tokens.Length() == 1)
               error("Unnamed marker in allele frequency file");
            if (info != NULL)
               need_blank_line |= info->AdjustFrequencies();
            info = GetMarkerInfo(tokens[1]);
            info->freq.Clear();
            info->freq.Push(0.0);
            next_allele = 1;
            break;
         case 'F' :
            if (info != NULL)
               for ( i = 1; i < tokens.Length(); i++)
                  {
                  buffer = next_allele++;

                  int allele = LoadAllele(info, buffer);

                  if (allele >= info->freq.Length())
                     {
                     old_max = info->freq.Length();
                     info->freq.Dimension(allele + 1);
                     for (j = old_max; j < allele; j++)
                        info->freq[j] = 0.0;
                     }

                  info->freq[allele] = tokens[i].AsDouble();
                  }
            break;
         case 'A' :
            if (info == NULL) continue;

            if (tokens.Length() != 3)
               error("Error reading named allele frequencies for locus %s\n"
                     "Lines with named alleles should have the format\n"
                     "   A allele_label allele_frequency\n\n"
                     "But the following line was read:\n%s\n",
                     (const char *) info->name, (const char *) buffer);

            allele_size = LoadAllele(info, tokens[1]);
            next_allele = atoi(tokens[1]) + 1;

            if (allele_size < 1)
               error("Error reading named allele frequencies for locus %s\n"
                     "An invalid allele label was encountered\n",
                     (const char *) info->name);

            if (allele_size >= info->freq.Length())
               {
               old_max = info->freq.Length();
               info->freq.Dimension(allele_size + 1);
               for (i = old_max; i < allele_size; i++)
                  info->freq[i] = 0.0;
               }

            info->freq[allele_size] = tokens[2];
            break;
         case 'E' :
            done = 1;
            break;
         default :
            error ("Problem in allele frequency file.\n"
                   "Lines in this file should be of two types:\n"
                   "  -- Marker name lines begin with an M\n"
                   "  -- Frequency lines begin with an F\n\n"
                   "However the following line is different:\n%s\n",
                   (const char *) buffer);
         }
      }

   if (info != NULL)
       need_blank_line |= info->AdjustFrequencies();

   if (need_blank_line) printf("\n");
   }

void PedigreeGlobals::LoadMarkerMap(const char * filename, bool filter)
   {
   IFILE f = ifopen(filename, "rb");
   if (f == NULL) return;
   LoadMarkerMap(f, filter);
   ifclose(f);
   }

void PedigreeGlobals::LoadMarkerMap(IFILE & input, bool filter)
   {
   String      buffer;
   StringArray tokens;
   bool first_pass = true;

   while (!ifeof(input))
      {
      buffer.ReadLine(input);

      tokens.Clear();
      tokens.AddTokens(buffer, WHITESPACE);

      if (tokens.Length() < 1) continue;

      if (first_pass)
         {
         sexSpecificMap = (tokens.Length() == 5);

         // if (sexSpecificMap)
         //   printf("\n  Found sex-specific map ...\n\n");

         first_pass = false;
         }

      if (tokens.Length() != 3 && !sexSpecificMap)
         error("Error reading map file\n"
               "Each line in this file should include 3 fields:\n"
               "CHROMOSOME, MARKER_NAME, and POSITION\n"
               "However the following line has %d fields\n%s\n",
               tokens.Length(), (const char *) buffer);

      if (tokens.Length() != 5 && sexSpecificMap)
         error("Error reading map file\n"
               "Each line in this file should include 5 fields:\n\n"
               "CHROMOSOME, MARKER_NAME, SEX_AVERAGED_POS, FEMALE_POS AND MALE_POS\n\n"
               "However the following line has %d fields\n%s\n",
               tokens.Length(), (const char *) buffer);

      bool previous_state = String::caseSensitive;
      String::caseSensitive = false;

      if ((tokens[0] == "CHR" || tokens[0] == "CHROMOSOME") &&
          (tokens[1] == "MARKER" || tokens[1] == "MARKER_NAME" || tokens[1] == "MRK") &&
          (tokens[2] == "KOSAMBI" || tokens[2] == "POS" || tokens[2] == "POSITION" ||
           tokens[2] == "SEX_AVERAGED_POS" || tokens[2] == "CM" || tokens[2] == "HALDANE"))
          continue;

      String::caseSensitive = previous_state;

      if (filter)
         if (LookupMarker(tokens[1]) < 0)
            continue;

      MarkerInfo * info = GetMarkerInfo(tokens[1]);

      int chr = (tokens[0][0] == 'x' || tokens[0][0] == 'X') ? 999 : (int) tokens[0];

      info->chromosome = chr;
      info->position = (double) tokens[2] * 0.01;

      if (sexSpecificMap)
         {
         char * flag;

         double female = strtod(tokens[3], &flag);
         if (*flag)
            error("In the map file, the female cM position for marker\n"
                  "%s is %s. This is not a valid number.",
                  (const char *) tokens[1], (const char *) tokens[3]);

         double male = strtod(tokens[4], &flag);
         if (*flag)
            error("In the map file, the male cM position for marker\n"
                  "%s is %s. This is not a valid number.",
                  (const char *) tokens[1], (const char *) tokens[4]);

         info->positionFemale = (double) female * 0.01;
         info->positionMale = (double) male * 0.01;
         }
      else
         info->positionFemale = info->positionMale = info->position;
      }

      if (sexSpecificMap) VerifySexSpecificOrder();
  }

void PedigreeGlobals::LoadBasepairMap(const char * filename)
   {
   IFILE f = ifopen(filename, "rb");
   if (f == NULL)
      error("The map file [%s] could not be opened\n\n"
            "Please check that the filename is correct and that the file is\n"
            "not being used by another program", filename);
   LoadBasepairMap(f);
   ifclose(f);
   }

void PedigreeGlobals::LoadBasepairMap(IFILE & input)
   {
   String      buffer;
   StringArray tokens;

   sexSpecificMap = false;

   while (!ifeof(input))
      {
      buffer.ReadLine(input);

      tokens.Clear();
      tokens.AddTokens(buffer, WHITESPACE);

      if (tokens.Length() < 1) continue;

      if (tokens.Length() != 3)
         error("Error reading map file\n"
               "Each line in this file should include 3 fields:\n"
               "CHROMOSOME, MARKER_NAME, and POSITION\n"
               "However the following line has %d fields\n%s\n",
               tokens.Length(), (const char *) buffer);

      bool previous_state = String::caseSensitive;
      String::caseSensitive = false;

      if ((tokens[0] == "CHR" || tokens[0] == "CHROMOSOME") &&
          (tokens[1] == "MARKER" || tokens[1] == "MARKER_NAME" || tokens[1] == "MRK") &&
          (tokens[2] == "BASEPAIR" || tokens[2] == "POS" || tokens[2] == "POSITION"))
          continue;

      String::caseSensitive = previous_state;
      
      MarkerInfo * info = GetMarkerInfo(tokens[1]);

      int chr = (tokens[0][0] == 'x' || tokens[0][0] == 'X') ? 999 : (int) tokens[0];

      info->chromosome = chr;
      info->position = (double) tokens[2];
      }
   }

int PedigreeGlobals::instanceCount = 0;

PedigreeGlobals::~PedigreeGlobals()
   {
   if (--instanceCount == 0 && markerInfoSize)
      {
      for (int i = 0; i < markerInfoCount; i++)
         delete markerInfo[i];
      delete [] markerInfo;
      delete [] markerInfoByInteger;
      }
   }

void PedigreeGlobals::WriteMapFile(const char * filename)
   {
   if (!MarkerPositionsAvailable())
      return;

   FILE * output = fopen(filename, "wt");

   if (output == NULL)
      error("Creating map file \"%s\"", filename);

   WriteMapFile(output);
   fclose(output);
   }

void PedigreeGlobals::WriteMapFile(FILE * output)
   {
   if (!sexSpecificMap)
      fprintf(output, "CHR  MARKER    POS\n");
   else
      fprintf(output, "CHR  MARKER    POS   POSF   POSM\n");

   for (int i = 0; i < markerInfoCount; i++)
      if (markerInfo[i]->chromosome != -1)
         if (!sexSpecificMap)
            fprintf(output, "%3d  %-10s  %g\n",
               markerInfo[i]->chromosome,
               (const char *) markerInfo[i]->name,
               markerInfo[i]->position * 100.0);
         else
            fprintf(output, "%3d %-10s %g  %g  %g\n",
               markerInfo[i]->chromosome,
               (const char *) markerInfo[i]->name,
               markerInfo[i]->position * 100.0,
               markerInfo[i]->positionFemale * 100.0,
               markerInfo[i]->positionMale * 100.0);
   }

void PedigreeGlobals::WriteFreqFile(const char * filename, bool old_format)
   {
   FILE * output = fopen(filename, "wt");

   if (output == NULL)
      error("Creating allele frequency file \"%s\"", filename);

   WriteFreqFile(output, old_format);
   fclose(output);
   }

void PedigreeGlobals::WriteFreqFile(FILE * output, bool old_format)
   {
   for (int i = 0; i < markerInfoCount; i++)
      {
      MarkerInfo * info = markerInfo[i];

      if (info->freq.Length() == 0) continue;

      fprintf(output, "M %s\n", (const char *) info->name);

      if (old_format && info->alleleLabels.Length() == 0)
         for (int j = 1; j < info->freq.Length(); j++)
            fprintf(output, "%s%.5f%s",
               j % 7 == 1 ? "F " : "", info->freq[j],
               j == info->freq.Length() - 1 ? "\n" : j % 7 == 0 ? "\n" : " ");
      else
         for (int j = 1; j < info->freq.Length(); j++)
            if (info->freq[j] > 1e-7)
               fprintf(output, "A %5s %.5f\n",
                      (const char *) info->GetAlleleLabel(j), info->freq[j]);
      }
   }

bool PedigreeGlobals::MarkerPositionsAvailable()
   {
   for (int i = 0; i < markerInfoCount; i++)
      if (markerInfo[i]->chromosome != -1)
         return true;

   return false;
   }

bool PedigreeGlobals::AlleleFrequenciesAvailable()
   {
   for (int i = 0; i < markerInfoCount; i++)
      if (markerInfo[i]->freq.Length() > 1)
         return true;

   return false;
   }

int PedigreeGlobals::LoadAllele(int marker, String & token)
   {
   return LoadAllele(GetMarkerInfo(marker), token);
   }

int PedigreeGlobals::LoadAllele(MarkerInfo * info, String & token)
   {
   int allele = info->GetAlleleNumber(token);

   if (allele >= 0) return allele;

   static unsigned char lookup[128];
   static bool init = false;

   if (!init)
      {
      init = true;

      for (int i = 0; i < 128; i++)
         lookup[i] = 0;

      for (int i = '1'; i <= '9'; i++)
         lookup[i] = 1;

      lookup[int('a')] = lookup[int('A')] = lookup[int('c')] = lookup[int('C')] = 2;
      lookup[int('g')] = lookup[int('G')] = lookup[int('t')] = lookup[int('T')] = 2;
      }

   int  first = token[0];
   bool goodstart = first > 0 && first < 128;

   if (token.Length() == 1 && goodstart && lookup[int(token[0])])
      return info->NewAllele(token);

   if (!goodstart || lookup[int(token[0])] != 1)
      return 0;

   int integer = token.AsInteger();
   token = integer;

   allele = info->GetAlleleNumber(token);

   if (allele > 0)
      return allele;

   if (integer <= 0) return 0;

   if (integer > 1000000)
      {
      static bool warn_user = true;

      if (warn_user)
         printf("Some allele numbers for marker %s are > 1000000\n"
                "All allele numbers >1000000 will be treated as missing\n\n",
                (const char *) info->name, warn_user = false);

      return 0;
      }

   return info->NewAllele(token);
   }


 
